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PREFACE 

The  work  reported  herein  was  conducted  by  the  Arnold  Engineering  Development 
Center  (AEDC),  Air  Force  Systems  Command  (AFSC),  under  Program  Element  65807F. 
The  results  were  obtained  by  ARO.  Inc..  AEDC  Division  (a  Sverdrup  Corporation 
Company),  operating  contractor  for  the  AEDC,  AFSC,  Arnold  Air  Force  Station. 
Tennessee.  The  work  was  done  under  ARO  Project  No.  E32-POA.  Elton  R.  Thompson 
was  the  Air  Force  project  manager.  The  manuscript  was  submitted  for  publication  on 
April  4,  1978. 
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1.0  INTRODUCTION 

Subsonic,  turbulent,  recirculating  flow  fields  exist  in  ducts  where  sudden  expansions 
or  contractions  occur,  behind  protuberances  such  as  codling  fins  and  vortex  generators,  in 
the  diverging  sections  of  wide  angled  diffusers,  within  combustors  and  cavities,  and.  to 
some  extent,  in  plenum  chambers  and  test  cells  of  propulsion  engine  test  facilities. 
Development  of  the  capability  to  compute  this  class  of  flows  will  greatly  enhance  the 
understanding  of  and  design  capability  for  recirculating  flow  configurations.  An  important 
example  of  this  class  of  flows  occurs  in  the  sudden  expansion  or  dump  combustors  that 
are  used  in  aircraft  and  missile  propulsion  systems.  A typical  flow  configuration  is  shown 
in  Fig.  1 which  illustrates  a recirculating  flow  created  by  coaxial,  two-stream  jet  mixing. 
The  operation  of  such  combustors  depends  on  the  stable,  complete  burning  that  occurs  in 
the  recirculation  zones  and  eddies  that  are  created  by  the  aerodynamic  forces  induced  by 
turbulent  mixing.  Hence,  it  is  important  that  the  flow  fields  within  such  combustors  be 
calculated  so  that  efficient,  clean  operation  can  be  obtained. 

Unfortunately,  recirculating  flows  are  mathematically  described  by  a coupled  set  of 
non-linear  or  quasi-linear  elliptic  partial  differential  equations.  The  set  can  contain  cither 
the  full  Navier-Stokes  equations,  or,  can  be  based  on  a pair  of  equations  describing  the 
distribution  of  stream  function  and  vorticity.  Solution  methods  for  such  equations  have 
been  developed  during  the  last  two  decades  based  on  finite  difference  techniques. 
However,  few,  if  any,  accurate  solutions  of  the  flow  field  shown  in  Fig.  1 have  been 
calculated  for  two  primary  reasons.  First,  solution  techniques  for  elliptic  flows  must  be 
based  on  iterative  procedures.  As  a result,  the  calculations  must  be  made  with  rather 
crude  finite  difference  meshes  to  keep  the  calculation  times  and  computer  storage 
requirements  within  practical  limits.  Second,  turbulent  viscosity  models  for  such  flows 
were  not  formulated  until  recently  and  still  remain  a topic  of  research.  Other, 
code-specific  problems  have  been  encountered.  For  example,  Schulz  (Refs.  1 and  2)  using 
the  elliptic  recirculating  flow  solution  method  provided  by  Gosman,  et  al.  (Ref.  3),  found 
that  predictions  of  recirculating  combustor  flows  of  acceptable  appearance  could  only  be 
obtained  after  extensive  trial  and  error  manipulation  of  three  of  the  principal  variables  of 
the  solution  procedure:  the  finite  difference  mesh  configuration  used  to  represent  the 
combustor  geometry,  the  boundary  condition  for  vorticity  at  the  lip  of  the  fuel  or 
oxidizer  nozzle,  and  the  constants  in  the  turbulent  viscosity  model.  Three  recirculating 
flows  of  the  type  illustrated  in  Fig.  1 were  calculated,  and  each  required  a somewhat 
different  combination  of  the  three  factors  for  the  best  results.  Another  significant  but 
unresolved  problem  identified  in  the  studies  reported  in  Refs.  1 and  2 was  that  the 
governing  equation  for  vorticity,  solved  by  the  code  provided  in  Ref.  3,  was  incomplete. 
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Terms  involving  the  derivatives  of  viscosity  arise  in  the  vorticity  transport  equation  when 
it  is  derived  from  the  Navier-Stokes  equations.  Gosman,  et  al.,  deleted  those  terms  from 
the  governing  equation  by  assuming  that  they  were  negligible  compared  with  all  the  other 
terms.  However,  the  effects  of  these  terms  on  the  calculated  flow  fields  have  not  been 
determined,  and  it  is  probable  that  they  are  not  negligible. 

Because  accurate  predictions  of  even  the  geometrically  simple  combustor  flows,  such 
as  Fig.  1,  have  not  been  made,  a study  was  undertaken  to  develop  a computer  program 
that  can  make  accurate  and  practical  calculations.  In  addition,  the  code  was  to  be  based 
on  the  complete  vorticity  transport  equation  for  turbulent  flows.  To  test  the  accuracy  of 
the  code,  four  turbulent  experiments  were  selected  to  compare  with  the  theory-an 
incompressible,  fully  developed,  turbulent  pipe  flow  and  three  cases  of  ducted, 
two-stream  jet  mixing  with  recirculation:  one  constant  density,  one  nonreacting 
hydrogen-air,  and  one  reacting  hydrogen-air. 

2.0  THEORY 


2.1  METHOD  OF  ANALYSIS 

The  method  used' in  the  present  study  to  compute  field  properties  of  recirculating 
flows  is  the  method  developed  by  Chien  (Ref.  4),  which  is  based  on  solving  a stream 
function  and  vorticity  formulation  of  the  Navier-Stokes  equations  for  fully  elliptic  flows. 
Chein  demonstrated  the  utility  of  solving  the  governing  systems  of  equations  by  a finite 
difference  method  using  coordinate  stretching  to  provide  numerical  resolution  where 
required  and  special  "decay  functions"  to  provide  numerically  stable  solutions. 

In  this  investigation,  the  Navier-Stokes  equations  were  manipulated  to  form  a 
vorticity  transport  equation  valid  for  variable-viscosity,  variable-density  flows.  Turbulent 
transport  models  and  a chemical  reaction  model  were  selected.  The  problem  was 
simplified  somewhat  by  considering  that  the  flows  were  to  occur  at  constant  static 
pressure,  and  in  fact,  the  four  experimental  flows  selected  for  comparison  with  the 
theory  do  occur  at  nearly  constant  static  pressure. 

2.2  DERIVATION  OF  THE  VORTICITY  AND  STREAM  FUNCTION 

GOVERNING  EQUATIONS 


The  Navier-Stokes  equations  for  two-dimensional  steady  flow,  taken  from  Ref.  5, 

are: 

pu  + pv  in.  = -dP  + rVf2du  _ 2 V ■ 

d Y,  <?r  dZ  di  L \ d7,  3 


+ 


(1) 
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where  V is  the  gradient  operator,  v is  the  total  velocity  vector,  and  6 is  a parameter  such  that 
6 = 0 for  plane  two-dimensional  flows  and  5 = 1 for  two-dimensional  flows  with  axial 
symmetry.  To  derive  a vorticity  transport  equation  for  two-dimensional  flows,  Eq.  (2) 
was  differentiated  with  respect  to  Z,  Eq.  (1)  was  differentiated  with  respect  to  r,  and  the 
resulting  equations  were  subtracted  one  from  the  other  to  obtain  a single  equation  which 
is  explicity  independent  of  pressure  terms.  The  definition  of  vorticity  in  two-dimensional 
flows. 


a = <hL  - in 

dZ  dr 


(3) 


was  then  introduced,  and  the  single  equation  remaining  from  three  operations 
on  the  Navier-Stokes  equation  was  manipulated  into  an  equation  for  conservation  of 
vorticity, 


d2o  + d2a  _ ( 

' £1 

. _ 2 dP  N 

da 

- (El  - 2 9p 

_ s\  da 

dZ2  dr2  ' 

> P 

p d 7.  / 

dZ 

V p p dr 

r / dr 

+ ri  + 

i 

ii  + s | 

Zpv 

+ J_  - J.'S 

la 

Lp  dZ2 

P 

dr2 

' pr 

fir  dr  r2  / 

J 

+ 2 r*V  ( 

dv 

- du  \ + 

d2p 

du  _ 92p  dv  + 

dp  1 / du.  + dv. 

- 

p |_dZdr  \ 

<9r 

dz) 

dZ2 

dr  ^r2  dz 

dZ  dr  VdZ  dr 

r / 

- d ( du 

97  dZ  \d7 

j. 

dv  T v8\ 
dr  r / 

]* 

1 I’d  <u2  + V2  ) 

2p  L dZ 

dp  _ d (u2  + v2) 
dr  dr 

dpi 

dzj 

= 0 (4) 

In  Eq.  (4),  viscosity  is  assumed  to  be  a general  function  for  Z and  r.  In  order  to  apply 
Eq.  (4)  to  turbulent  flows,  the  only  assumption  required  is  that  the  viscosity  (p)  is  an 
"effective"  or  "eddy"  viscosity. 

The  stream  function  equation  was  derived  by  defining  a stream  function  (&)  such 

that 


d \ft  = pur^  dr  — pvr^  dZ 


(5) 


7 


AEDC-TR-78-29 


Then  the  definition  (Eq.  (5))  was  introduced  into  Eq.  (3),  and  the  result  was 

manipulated  to  form 

<?y  + - 8 ^ - I # - 1.  dp.  + pr^Q  = 0 (6) 

S'/2  dr 2 T dr  PS'/.  67.  P dr  dr 

Equation  (6)  is  the  governing  equation  for  the  stream  function.  It  is  easy  to  show  that 
the  stream  function  (i //)  automatically  and  identically  satisfied  the  mass  conservation 
equation  for  steady  flow.  i.e.. 


jL  (pu)  + JL  (pvr^)  =0  (j) 

d7  dr  K > 

2.3  TURBULENT  VISCOSITY  MODEL 


2.3.1  Constant  Viscosity 

In  this  model,  the  turbulent  viscosity  is  set  equal  to  a selected  constant, 

Pi  = Prkf  (8) 

The  use  of  Eq.  (8)  effectively  treats  the  turbulent  flow  as  a laminar  flow. 

2.3.2  Algebraically  Prescribed  Viscosity 


Schulz  (Ref.  1 and  2)  developed  the  following  model  to  correlate  theoretical 
predictions  of  recirculating  combustor-type  flows  with  experimental  data: 

puj  Rj  [1  + Z/(275llj)]  (9) 

where  is  a constant  whose  value  is  0.0285.  Equation  (9)  was  developed  for 
two-stream,  coaxial  jet  mixing  with  recirculation  when  the  duct  to  jet  radius  ratio  is 
greater  than  ten.  There  is  no  physical  justification  for  this  model  except  that  fully 
developed,  turbulent,  free  jets  can  be  shown  to  require  a constant  viscosity  of  magnitude 
(k^piijRj).  Numerical  experimentation  by  Schulz  showed  that,  for  the  recirculating  flows 
of  this  type,  a linearly  increasing  effective  viscosity  was  required  for  the  best  prediction. 

2.3.3  Prandtl  Mixing  Length  Viscosity  Model 


The  mixing  length  model  developed  by  Prandtl  is  well  described  by  Schlichting  (Ref. 
6).  It  has  been  used  to  calculate  turbulent  jets,  wakes,  boundary  layers  and  pipe  flows.  In 
Prandtl's  model, 


Pi  = Pi2 


(10) 
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where  2 is  a length  scale  called  the  Prandtl  mixing  length.  For  fully  developed  turbulent 
flow  in  smooth  walled  pipes,  Schlichting  gives 

i = [o.l4  - 0.08  (r/Ru)2  - 0.06  Cr/ R D) 4 ] RD  (11) 


2.3.4  Turbulent  Kinetic  Energy-Dissipation  Viscosity  Model 


This  model  for  turbulent  viscosity  was  developed  by  Spalding  and  his  co-workers 
and  is  described  by  Launder  and  Spalding  (Ref.  7).  It  has  been  applied  to  several  types  of 
turbulent  flows,  including  flows  with  recirculation  (Ref.  8)  but  with  only  partial  success. 
In  this  viscosity  formulation,  the  turbulent  viscosity  is  given  as 

pt-cMpKa/e  (12) 

where  is  a constant  with  the  value  0.09.  To  determine  the  spatial  distribution  of  K 
and  e in  high  turbulent  Reynolds  number  flows,  two  coupled  partial  differential  transport 
equations  must  be  solved: 


(pu  K)  -t  (p  v K) 

dZ  r dr 

-k 

(JL 

v8  dK\ 

dz  ) 

- jL 

dr 

(JL  r8  dK'S  - rSS„  = 0 
W dr  ) K 

(13) 

JL  (p  ur^f)  +_d.(pvr^e) 

SZ  dr 

d 

- dz 

(<L 

\ 

3 K) 

- JL 

dr 

(JLt8  \ S = 0 

\*e  dr  ) f 

(14) 

where  the  parameters  Or 

and 

o€ 

are  turbulent 

Schmidt  numbers  for  K 

and  e 

respectively,  having  the  values  pr  = 1.0  and  ae  = 1.3. 


Sk  and  Se  are  source  terms  for  K and  e given  by 


where  C£l  and  Cg2  are  constants  having  the  values  of  1 .45  and  2.0,  respectively. 


The  constant  effective  viscosity  model  was  used  in  the  evaluation  and  checkout  stage 
of  the  program  development.  The  Prandtl  mixing  length  model  and  the  K-e  model  were 
used  to  calculate  the  turbulent  pipe  flow,  while  the  algebraic  viscosity  model  was  used  in 
the  calculation  of  the  jet  mixing  with  recirculation  flows. 
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2.4  EQUAf  ION  OF  STATE 

For  incompressible  flows,  the  density  was  set  to  an  appropriate  constant  value,  but 
for  the  nonreactive  and  reactive  calculations,  the  density  at  each  point  in  the  flow  was 
calculated  from  the  ideal-gas  equation  of  state 

r _ p » (17) 


For  the  nonreactive  calculation,  both  pressure  and  temperature  were  constant,  but  the 
molecular  weight  (W)  had  to  be  calculated  at  each  point  in  the  flow.  The  molecular 
weight  was  computed  from  the  distributions  of  the  mass  fractions  of  hydrogen  and  air. 
Once  the  local  value  of  the  mass  fraction  of  air  (C)  was  known,  the  molecular  weight  was 
determined  from 


W = 1 //__£_  i-  i - c\ 

f\  28.96  2 ) (18) 

To  determine  the  spatial  distribution  of  C,  a model  transport  equation,  valid  for 
turbulent  flows,  was  solved  simultaneously  with  the  governing  equations.  The  equation, 
obtained  from  Gosman,  et  al.  (Ref.  3),  is 


JL  (pur5C)  - JL  (pv rSC)  - JL  ( JL 

d'/.  dr  dZxO'c 


rS  K\  - JL  ( JL  rSdC\  = 0 

dZ  ) dr  \<7c  dr  J 


(19) 


For  the  reacting  flow  calculations,  the  pressure  was  assumed  constant;  however,  both 
the  temperature  and  the  molecular  weight  had  to  be  determined  at  each  point  of  the 
flow  field.  To  determine  their- variations,  a simplified  chemical  reaction  model  and  a thin 
flame  sheet  model  were  used.  This  approach  was  suggested  by  Gosman.  ct  al.  (Ref.  3).  In 
the  present  study,  the  following  chemical  reaction  model  for  the  hydrogen-oxygen- 
nitrogen  system  was  adopted : 

l Kg  Fuel  + i Kg  Oxidizer  + jKg  Inert  -«  (i  + 1)  Kg  Products  + jKg  Inert  where  for 

the  h\drogen-air  reaction,  i = 8 and  j = 26.3.  (20) 

The  thin  flame  sheet  model  was  developed  to  calculate  the  reaction  zone  of  initially 
separate  fuel  and  oxidizer  streams.  It  is  assumed  that  the  fuel  and  oxidizer  streams  are 
separated  by  a thin  flame  sheet  such  that,  on  the  fuel  side  of  the  flame  sheet,  only  the 
mass  fractions  of  fuel,  products,  and  inert  species  exist;  whereas,  on  the  oxidizer  side  of 
the  flame  sheet,  only  the  mass  fractions  of  oxidizer,  products,  and  inert  species  exist.  The 
location  of  the  thin  flame  sheet  is  the  locus  of  points  where  the  mass  fractions  of  both 
the  fuel  and  oxidizer  have  disappeared.  To  predict  the  spatial  variations  of  the  specie 
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mass  fractions,  a transport  equation  for  the  chemical  reaction  parameter  C was  derived  by 
following  the  procedure  outlined  by  Gosman,  et  al.  (Ref.  3).  The  equation  is  the  same  as 
the  mass  fraction  transport  equation  for  the  nonreactive  flow: 


jL  (p  u C)  + 

37  dr 


(pv rSC)  - A (JL  Ts  d_C\  _ A (JL  rSd£\ 
d7  \CTc  3Z/  dr  V^c  dr) 


0 


(21) 


but  for  the  reactive  flow 


C = (C-  Cfi)/(Q 


- 


(22) 


where  FI  and  01  are  the  fuel  and  oxidizer  inlet  conditions,  respectively,  and 

c = mqj  - imF  (23)' 

The  equation  for  C is  solved  simultaneously  with  the  governing  equations.  Once  C had 
been  found,  locally,  the  mass  fractions  of  the  species  were  known.  Specifically  for  the 
hydrogen-air  reaction  on  the  fuel  side  of  the  flame  sheet 

m,  = 0.7678  C 
mcp  = 0 

mF  = 1 - 1.029025  C 

m = 0.261225  C (24) 

while  on  the  oxidizer  side  of  the  flame  sheet 


m,  = 0.7678  C 

= 8.2322  C - 8 

m F = 0 .0 

mp  = 9(1  - C)  (25) 

Thus,  the  local  value  of  the  molecular  weight  can  be  determined  indirectly  as  a 
function  of  C,  since 


W 


2 


!*P  + 

32 


(26) 
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To  determine  the  temperature  field,  the  following  assumptions  were  made: 

1.  The  kinetic  energy  of  the  mean  flow  and  the  turbulence  are  small 
compared  to  the  chemical  energy  of  the  hydrogen-air  reaction. 

2.  The  specific  heats  of  the  chemical  species  are  constant. 

3.  The  total  enthalpy  of  the  flow  which  includes  both  the  sensible  and 
chemical  energy  is  everywhere  constant.  Also  the  walls  of  the  duct  are 
adiabatic. 


By  making  these  assumptions,  the  static  temperature  at  each  point  in  the  flow  field  is 
given  by 


T = T 


R F K 


+ m. 


11  n ff  7 <CpK 


m. 


CPlmI 


C m 
P<* 


+ c 


Pp  P 


(27) 


The  reference  temperature  was  assumed  to  be  the  static  temperature  of  the  inlet 
hydrogen  and  air  streams  which  were  equal.  The  reference  enthalpy  was  selected  such 
that  the  maximum  temperature  in  the  flame  sheet  was  the  stoichiometric  temperature  for 
the  hydrogen-air  reaction.  The  values  of  these  parameters  were  specified  as 

Tref  = 294  °K 

IIRKF  = 12317656.9  J/Kg 

C.  = 16099.0  J/Kg-°K 
ff 

Cpi  = 1173.9  J/Kg-°K 

Cn  = 1048.1  J/Kg-flt 

?<f> 

Cpp  = 2431.6  J/Kg-°K  (28) 

For  these  conditions,  the  stoichiometric  flame  temperature  for  hydrogen  and  air  is 


T STOIC  = 2388‘9 


(29) 


Thus,  the  density  of  the  reacting  hydrogen-air  flow  was  determined  (indirectly  as  a 
function  of  C)  according  to  the  equation 


P 


_JP  ^ ( 294  + , 12337656.9  mp  

32  7 V (1  173.9  nij  + 2431. 6 mp  + 16099mF  + 1048.1  m ) 


(30) 
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2.5  BOUNDARY  CONDITIONS 

The  boundary  conditions  for  the  dependent  variables  were  different  for  the 
turbulent  pipe  flow  and  the  jet-mixing  flows.  The  boundary  conditions  are  specified  below 
for  both  classes  of  turbulent  flow. 

2.5.1  Turbulent  Pipe  Flow 

2.5.1. 1 Vorticity 

Across  the  inlet  of  the  pipe,  the  velocity  was  assumed  uniform  so  that  12  = 0.  Along 
the  wall,  or  rather  at  the  first  mesh  point  off  the  wall,  a slightly  modified  form  of  a 
formula  used  by  Chien  (Ref.  4)  was  adopted, 

°I.JN  ~ -Ol.JK-2  + - ('''J"  L"''1’1-2)  <31> 

Along  the  axis  of  the  pipe,  12  = 0,  while  across  the  exit  of  the  pipe,  the  condition  312/3 Z 
= 0 was  imposed. 

2.5.1. 2 Stream  Function 

For  uniform  inlet  flow,  the  stream  function  distribution  is  i£(r)  = 1/2  pur2.  Along 
the  wall  and  axis  of  the  pipe,  the  stream  function  is  constant.  At  the  outlet  of  the  pipe, 
the  condition  3i///3Z  = 0 was  imposed. 

2.5. 1.3  Turbulent  Kinetic  Energy  and  Dissipation 

Across  the  inlet  of  the  pipe,  K and  e were  set  to  a small  constant  (10‘6).  Along  the 
wall,  the  formulas  developed  by  Chien  (Ref.  4)  were  applied,  i.e., 

kjk  - (32) 

* JN  " PjK  V"  %J*  (33) 

where  V*  is  the  local  value  of  the  friction  velocity  defined  asyVw/p.  In  some  calculations,  V* 
was  assigned,  based  on  experimental  data;  whereas  in  other  computations,  V*  was 
calculated  using  the  iterative  law-of-the-wall  matching  procedure  of  Chien  (Ref.  4).  Across 
the  pipe  outlet,  the  condition  3K/3Z  = 3e/3Z  = 0 was  imposed;  while  along  the  axis  of 
the  pipe,  3K/3r  = 3e/3r  = 0. 
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2.5.2  Ducted  Recirculating  Flow 

2.5.2. 1 Vorticity 

Across  the  jet  inlet,  from  the  axis  to  the  nozzle  lip,  the  vorticity  distribution  was 
assigned  by  assuming  that  Q.  = -Au/Ar.  The  necessary  values  for  Au/Ar  were  obtained 
from  a 15-point  central  difference  analysis  of  the  velocity  profile  shown  in  Fig.  2.  This 
profile  was  obtained  from  the  smoothed  experimental  jet  profile  of  Ref.  2.  The  velocity, 
vorticity,  and  stream  function  distributions  were  truncated  at  the  radial  location  r = Rj  - 
5*  where  5*  is  the  displacement  thickness  of  the  measured  velocity  profile.  In  principle, 
the  truncated  part  of  velocity  and  velocity  profile  should  have  a negligible  effect  on  the 
downstream  flow  field,  and  therefore,  it  was  assumed  that  it  was  unnecessary  to 
accurately  model  this  part  of  the  jet  exit  flow  field. 

Across  the  remainder  of  the  inlet  plane  from  jet  nozzle  lip  to  the  duct  wall,  along 
the  wall,  and  along  the  centerline,  the  vorticity  was  set  equal  to  zero.  Across  the  duct 
outlet  plane  the  condition  d£2/dZ  = 0 was  imposed. 

2.5.2. 2 Stream  Function 

Across  the  inlet  of  the  duct,  from  the  axis  to  the  jet  nozzle  lip,  the  stream  function 
distribution  was  obtained  by  assuming  constant  jet  flow  density  and  numerically 
integrating  the  velocity  profile  (Fig.  2).  The  integration  was  carried  out  such  that  the 
relation, 


u 


Pr 


Ar 


(34) 


would  reproduce  the  velocity  profile  from  the  stream  function  profile.  The 
self-consistency  of  the  velocity  and  stream-function  profiles  are  indicated  in  Fig.  2 by  the 
excellent  agreement  of  the  calculated  velocity  with  the  given  profile. 


The  stream-function  distribution  over  the  remainder  of  the  duct  inlet  from  the  jet 
nozzle  lip  to  the  wall  was  determined  by  assuming  constant  secondary  stream  density  and 
velocity,  so  that 


J P s us^2  + ^JET  LIP  : r > Rj  (35) 

Along  the  wall  and  along  the  axis  of  the  duct,  the  stream-function  is  constant.  Across  the 
outlet  of  the  duct,  the  condition  b\}jjbZ  = 0 was  imposed. 
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2. 5. 2. 3 Species  Parameter 

Across  the  inlet  of  the  duct,  from  the  axis  to  the  jet  nozzle  lip,  C = 1.0,  and  from 
the  nozzle  lip  to  the  wall,  C = 0.0.  Along  the  wall  and  along  the  axis,  the  condition 
3C/9r  = 0 was  imposed,  while  across  the  outlet  plane  the  condition  3C/3Z  = 0 was 
imposed. 


3.0  NUMERICAL  SOLUTION  PROCEDURE 

The  numerical  solution  procedure  used  in  the  present  study  was  developed  and 
reported  by  Chien  (Refs.  4,  9,  and  10).  The  numerical  solution  procedure  is  written  in 
terms  of  stretched  coordinates.  Chien  used  coordinate  stretching  to  transform  a physical 
domain  enclosed  by  diverging  diffuser  walls  into  a rectangular  computational  domain  and 
to  provide  numerical  resolution  in  the  boundary  layer.  However,  in  the  present  study, 
coordinate  stretching  is  used  to  stretch  the  flow  field  near  the  central  inlet  in  order  to 
improve  the  numerical  resolution  of  the  solution.  The  coordinate  stretching  required  in 
the  present  study  was  provided  by  simple  transformation  functions  of  the  form 

e = B(z) 

r)  = f(r)  (36) 

This  family  includes,  for  example,  power  law  and  logarithmic  stretching  functions.  The 
transformation  factors  are 


8_  = H d_ 

87.  87  8£ 

8_  = JL 

dr  dr  drj 

JL  = ^_J.  + (*£\2  JL 

87 2 <9/ 2 M \dz  J 

JL  = dLn  jL  - (Ll\2  JL 

dr 2 dr2  dr}  \8  r)  dr} 2 

82  _ <%_  LL  d2 

87.8 r 87  dr  dtjdr} 


(37) 


The  coordinate  transformations  (Eqs.  (36)  and  (37))  were  applied  to  the  set  of  governing 
partial  differential  equations,  and  then  each  equation  of  the  set  was  rearranged  to  the 
"standard  form"  proposed  by  Chien  (Ref.  4) 
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al 


- a2  . ^ - bL  , ^ - b2  , ^ + d , = 0 

de  9 9 ^ dv  * 


(38) 


The  coefficients  (al^,  a2$.  bl^  and  b2^)  and  the  source  term  (dp)  are  listed  in  Appendix 
A for  each  of  the  parameters  0 = S2,  i p,  C,  K,  and  e. 


In  the  "decay  function  formulation,"  functions  are  introduced  into  the  finite 
•difference  equations  that  correspond  to  the  governing  equations  to  make  the  Finite 
difference  equations  exact.  The  decay  functions  are  obtained  from  local  exact  solutions 
of  the  governing  differential  equations.  If  the  governing  equations  are  partial  differential 
equations  or  if  the  coefficient  al^.  a2^,  bl<£,  or  b2p  vary  as  functions  of  coordinate 
location,  then  the  decay  functions  provided  by  Chien's  theory  are  themselves 
approximate.  Interestingly,  although  all  finite  difference  equations  written  in  analogy  to 
Eq.  (38)  are  cast  in  central  or  space  centered  finite  differences,  the  decay  function 
formulation  approaches  the  classical  upwind  one-sided  difference  formulation  when  the 
local  cell  Reynolds  number  becomes  much  greater  than  two. 


The  local  cell  Reynolds  numbers  are  defined  as 


R*f  - “d  -v 


and 


' = b2^A,?/a20 
The  decay  functions,  formulated  in  the  approximate  form  presented  by  Chien,  are 

G0£  = 1 - 0.0625  |<  2 


(39) 


h.  n -> 


>2 


(40) 


and 


G*  = 1 - 0.0625  Rj 


R*  1 <2 


'’I 


J 


> 2 
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Each  governing  equation  is  written  in  the  space-centered  finite  difference  form 

Hjl  . c. 


’*€ 


1 

n 

i 

+ 

i— » 

«— i 

+ 

■o" 

20ij\ 

\e 

) 

- aj±  [y<v.T-i+<vi-i-2<?u>\  - / % \ \ 

L\  V ) \ ^ J \ 2^  )_ 


+ \'j  = ° (4D 

This  equation  is  solved  algebraically  for  0i(j  which  results  in  the  equation 

^i, J = c0,+  1 0i  + i, J + 0i-i. J + C</>J+1  0i,j  + i + 0i.j-i  + D0,  j 

(42) 

where  the  coefficients  are 


bV2A<l/s 

V.  - K/4f2G*{  - kv“f]/% 

C*Jtl  - 

<V,  ■ + bv24*]/S 

CV  = 2[al0M^G0£  * aV/A,? 2 ] (43) 

and  the  source  term  is 

% ■ (+*> 

Equation  (42)  was  used  as  a successive  substitution  formula  for  An  iterative  field 
point  subroutine  was  developed  to  obtain  convergent  solution  for  each  variable  (J2,  ijj,  C, 
K,  and  e)  in  a point-by-point  manner.  The  rate  of  convergence  was  controlled  by  an 
appropriate  relaxation  factor. 
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4.0  COMPARISON  OF  THEORY  AND  EXPERIMENT 
4.1  TURBULENT  PIPE  FLOW 

The  numerical  algorithm  was  checked  by  using  the  constant  effective  viscosity  model 
to  predict  fully  developed  laminar  pipe  flow.  The  predictions  agreed  exactly  with  the 
known  analytic  solution  (Ref.  6).  Then  a series  of  calculations  was  made  of  a turbulent 
pipe  flow  investigated  and  reported  by  Laufer  (Ref.  11).  Distributions  of  axial  velocity 
and  turbulent  kinetic  energy  across  the  flow  in  a smooth-walled  pipe  were  derived  from 
Laufer's  data  for  a Reynolds  number  of  500.000.  which  served  as  the  basis  for 
comparison  with  the  theory. 

The  first  calculation  was  made  by  using  the  experimental  data  to  establish  the 
boundary  conditions  for  velocity  and  friction  velocity.  The  boundary  conditions  were 
applied  at  the  radial  location  7 = 0.05  because  the  viscosity  models  are  not  valid  in  the 
laminar  sublayer  which  lies  next  to  the  wall.  The  Prandtl  mixing  length  viscosity  model 
was  applied  to  make  this  calculation.  Comparison  of  theory  and  experiment  is  shown  in 
Fig.  3.  The  velocity  field  is  excellently  predicted  but  the  turbulent  kinetic  energy 
distribution  near  the  wall  is  lower  than  experiment.  The  next  calculation  was  made  by- 
applying  a law-of-the-wall  matching  procedure  used  by  Chien  (Ref.  4)  for  assigning  values 
of  velocity  outside  the  laminar  sublayer.  In  this  procedure,  the  friction  velocity  and  the 
axial  velocity  at  the  boundary  point  (7  = 0.05)  are  obtained  during  the  course  of  the 
solution  by  using  the  velocity  at  the  second  point  from  the  wall  (r  = 0.105)  in  the 
law-of-the-wall  relationship  . a 


= 2.5  loge  (V*  y/v)  + 5.5 

to  solve  for  V*.  Then,  the  velocity  at  the  boundary  point  is  calculated  from  Eq.  45.  The 
results  of  this  calculation  are  shown  in  Fig.  4.  Again,  the  axial  velocity  distribution  is 
well  predicted,  but  the  turbulent  kinetic  energy  distribution  falls  below  the  experiment. 
Thus,  the  boundary  value  for  K is  too  small,  since  with  = 0.09 

K ,n  = V*  2 0.09  = 3.33V*2  (46) 

hence,  as  shown  in  Fig.  4,  0.2  K/V*2  = 0.67  at  the  matching  point,  7 = 0.05.  But  the 
experimental  data  imply  that  0.2  Kw/V*2  = 0.9  which  requires  that  C p = 0.05  at  the 
wall.  So,  as  Chien  (Ref.  4)  and  others  have  pointed  out,  C ^ must  vary  across  the 
turbulent  pipe  flow,  tending  to  about  0.05  at  the  boundary  point.  Therefore,  for  the 
following  calculations,  C ^ was  allowed  to  vary  linearly  from  0.09,  at  the  radial  location  r 
- 0.25,  to  0.05  at  the  wall,  7 = 0.  Both  the  Prandtl  mixing  length  viscosity  model  and  the 
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K-e  viscosity  model  were  used  to  calculate  the  flow  by  applying  the  law-of-the-wall 
matching  procedure.  The  results  are  shown  in  Fig.  5.  Both  models  provided  good 
predictions  of  the  axial  velocity  distributions,  although  the  Prandtl  mixing  length  model 
does  best,  and  both  models  predicted  about  the  same  turbulent  kinetic  energy 
distributions  which  are  in  good  agreement  with  the  experimental  data.  Also,  both  models 
predicted  the  wall  friction  velocity  to  within  three  percent. 

The  Reynolds  number  for  all  the  calculations  was  5.2  x 105  based  on  the  pipe 
diameter;  the  calculations  with  a 20  by  20  uniform  mesh  took  less  than  five  minutes 
using  an  IBM  370/165  computer. 

4.2  DUCTED,  TWO-STREAM,  COAXIAL  JET  MIXING  WITH  RECIRCULATION 


4.2.1  Constant  Density  Jet  Mixing 


An  experiment  reported  by  Barchilon  and  Curtet  (Ref.  12)  was  selected  as  a test 
case  for  this  type  flow.  The  data  reported  were  taken  in  an  experimental  configuration 
similar  to  that  shown  in  Fig.  1.  To  make  calculations  of  this  flow,  the  algebraic  viscosity 
model  (Eq.  (9))  was  applied  together  with  the  boundary  conditions  described  in  Section 
2.5.2.  It  should  be  noted  that  the  experimental  data  were  obtained  in  a w'ater  tunnel 
where  the  jet  velocity  was  about  12.8  m/sec,  and  the  duct-to-jet  radius  ratio  was  13.33. 
In  the  calculations,  the  fluid  is  air,  the  jet  velocity  is  205.7  m/sec,  and  the  duct-to-jet 
radius  ratio  is  10.  However,  both  flows  have  the  same  Craya-Curtct  number  and, 
therefore,  should  be  similar  according  to  theory  (Ref.  13). 


Figure  6 shows  a comparison  between  experimental  and  theoretical  distributions  of 
axial  velocity  on  the  centerline  of  the  duct.  While  the  theoretical  predictions  appear  to 
decay  at  the  proper  rate,  the  length  of  the  potential  core  of  the  central  jet  is  not  well 
predicted.  Moreover,  the  predicted  velocity  slightly  overshoots  the  inlet  or  boundary 
value  for  jet  velocity.  On  the  other  hand,  the  distributions  of  axial  velocity  across  the 
duct  are  in  reasonably  good  agreement  with  the  experiment,  as  shown  in  Fig.  7. 


As  was  stated  in  the  Introduction,  one  objective  of  the  present  investigation  was  to 
determine  whether  the  vorticity  source  terms  that  have  been  neglected  in  previous 
recirculating  flow  calculations  (e.g.,  Refs.  1,  2,  3,  and  8)  are  significant.  For  the  constant 
density  jet  calculations,  numerical  experiments  were  performed  by  repeating  the 
calculation- of  the  Barchilon  and  Curtet  experiment  but  without  the  terms 


sn  = 
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in  the  vorticity  source  term  in  Eq.  (4).  The  calculated  velocity  field  is  also  shown  in  Figs. 
6 and  7 where  the  results  are  nearly  indistinguishable.  Thus,  for  constant-density  flows 
with  slowly  varying  eddy  viscosity,  i.e..  for  nearly  laminar  incompressible  flows,  the 
vorticity  source  terms  (Eq.  (47))  may  be  neglected. 

These  calculations  required  about  five  minutes  using  a 20  by  50  mesh.  The  mesh 
was  uniform  in  the  radial  direction  and  slightly  stretched  in  the  axial  direction. 

4.2.2  Variable-Density,  Nonreactive  Jet  Mixing 

An  experiment  reported  by  Schulz  (Refs.  1 and  2)  of  ducted,  coaxial,  two-stream 
variable-density  jet  mixing  was  selected  as  a test  case.  The  flow  field  corresponds  to  that 
shown  in  Fig.  1.  The  variable-density  flow  was  isothermal,  and  the  density  variations  were 
due  to  the  molecular  weight  difference  between  the  central  and  annular  streams.  The 
mixing  in  this  experiment  occurred  between  a central  air  jet  and  a surrounding  stream  of 
hydrogen  in  which  the  air  jet  inlet  velocity  was  about  206  m/sec.  The  ratio  of  the  duct 
radius  to  the  jet  radius  was  10. 

The  algebraic  viscosity  model  (Eq.  (9))  was  again  used  together  with  the  same  mesh 
and  boundary  conditions  used  to  calculate  the  constant  density  case.  The  complete 
vorticity  source  term  was  included  in  the  vorticity  transport  equation.  Equation  (19)  was 
solved  for  the  mass  fraction  of  air.  Figure  8 shows  a comparison  of  experimental  and 
theoretical  distributions  of  hydrogen  mass  fraction  along  the  centerline  and  the  wall  of 
the  duct  and  the  distribution  of  axial  velocity  along  the  centerline  of  the  duct.  The 
predicted  velocity  distribution  is  in  satisfactory  agreement  with  the  experiment.  But,  the 
predicted  hydrogen  mass  fractions  show  greater  diffusion  and  spreading  than  experiment. 
Also  the  predicted  fully  mixed  mass  fraction  is  eight  percent  higher  than  experiment.  On 
the  other  hand,  the  predicted  axial  velocity  distributions  across  the  duct  are  in  good 
agreement  with  the  data  (Fig.  9).  The  hydrogen  mass  fraction  distributions  at  the 
corresponding  axial  positions  are  compared  with  data  in  Fig.  10,  which  also  shows  that 
the  predicted  hydrogen  mass  fraction  field  has  greater  spreading  and  diffusion  than  the 
experiment  indicates.  Thus,  it  would  appear  that,  for  the  variable-density  flow,  a better 
turbulent  viscosity  model  is  required. 

The  effect  of  deleting  the  viscosity-dependent  vorticity  source  terms  (Eq.  (47))  in 
the  variable-density  calculation  is  also  shown  in  Figs.  8 and  9.  Figure  8 shows  that  by 
deleting  the  terms  (Eq.  (47))  significant  velocity  overshoot  occurs  in  the  near  field.  Thus, 
for  this  specific  flow,  these  terms  have  a significant  effect  on  the  predicted  flow  field.  On 
the  other  hand,  comparison  of  the  axial  velocity  distributions  across  the  duct  shows  that 
the  terms  have  a relatively  unimportant  effect  on  the  shapes  of  the  far  field  velocity 
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profiles.  Therefore,  the  complete  vorticity  source  term  should  be  incorporated  in  variable 
density  calculations  because  of  their  effect  on  the  near  field. 

Other  numerical  experimentation  showed  that  the  inlet  plane  distribution  of 
viscosity  can  also  affect  the  calculated  flow  field  one  mesh  line  downstream  because  it 
enters  the  calculation  through  the  viscosity-dependent  terms  of  Eq.  (47).  Thus,  by 
including  the  full  vorticity  source  term  in  the  calculation  procedure,  all  of  the  inlet  plane 
boundary  values  must  be  consistently  described,  and  this,  in  a sense,  forces  the  solution 
procedure  to  be  physically  realistic.  Calculation  time  for  this  case  was  approximately  20 
min  on  the  IBM  370/165. 

4.2.3  Variable-Density,  Reacting  Jet  Mixing 

A real  combustor  flow  field  is  characterized  by  regions  wherein  large  heat  release 
rates  from  the  burning  gases  cause  both  molecular  weight  and  temperature  gradients.  So 
that,  in  turbulent  reacting  flow,  the  time-averaged  density  depends  on  the  time-averaged 
temperature  and  the  time-averaged  molecular  weight.  The  process  can  be  represented 
either  as 


p = (PW/RgT)  (48) 

or  as 

p = p vv  / R G T (49) 

In  the  present  study,  the  latter  definition  has  been  adopted,  as  explained  in  Section  2.4. 
The  approach  was  applied  to  calculate  the  reacting,  recirculating  combustor  flow  reported 
by  Chriss  (Ref.  14).  The  data  obtained  by  Chriss  were  taken  in  the  configuration  shown 
in  Fig.  1.  The  mixing  and  chemical  reaction  was  between  a central  air  stream  and  a 
surrounding  hydrogen  stream  having  the  same  velocity  and  radius  ratios  as  those  of  the 
previously  discussed  nonreactive  ease. 

Theoretical  predictions  were  made  without  changing  the  mesh,  the  inlet  or  boundary 
conditions,  or  the  viscosity  model  from  those  applied  in  the  nonreactive  flows.  The 
equation  for  the  species  parameter  (Eq.  (21))  was  added  to  the  set  of  governing 
equations.  The  theoretical  predictions  are  compared  with  experiment  in  Figs.  11,  12,  and 
13.  While  the  hydrogen  at  any  point  of  the  flow  may  be  contained  in  either  water  vapor 
or  atomic  or  molecular  hydrogen,  the  experimental  data  do  not  enable  the  distribution  of 
the  hydrogen  among  the  molecular  species  to  be  determined.  Thus,  the  hydrogen  mass 
fractions  indicated  in  the  figures  is  the  total  or  atomic  mass  fraction.  Also,  as  shown  in 
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Figs.  11  and  12,  there  are  two  velocity  curves  marked  "experimental."  One  velocity  curve 
represents  velocity  measured  by  laser  velocimetry.  The  other  "experimental"  velocity 
curves  were  calculated  from  the  measured  pitot  and  static  pressures  and  the  hydrogen 
atom  mass  fraction  by  assuming  that  chemical  equilibrium  existed  at  each  point  in  the 
flow  Field. 

The  predicted  velocity  distribution  along  the  centerline  of  the  duct  shows  a very 
large  overshoot,  which,  however,  falls  between  the  two  "experimental"  curves.  This 
makes  it  difficult  to  determine  whether  the  chemistry  or  the  turbulent  viscosity  model  is 
in  error.  Another,  more  significant  problem  is  that  the  theoretically  predicted  axial 
distributions  of  hydrogen  mass  fraction  exceeds  the  experiment  in  the  far  field  by  75 
percent.  Thus,  the  numerical  solutions  are  not  preserving  the  species  mass,  probably 
because  of  the  large  overshoot  of  the  velocity  Field.  In  order  to  resolve  the  problem, 
numerical  experimentation  would  be  required  beyond  the  scope  of  the  present  study.  The 
axial  velocity  profiles  across  the  flow  at  three  axial  stations  are  compared  with 
experiment  in  Fig.  12.  They  are  not  as  well  predicted  as  the  nonreactive  case  showing 
greater  spreading  than  the  experimental  data.  The  hydrogen  mass  fraction  profiles  across 
the  duct  arc  compared  with  experiment  in  Fig.  13.  The  shapes  of  the  calculated  hydrogen 
atom  fraction  profiles  are  fairly  well  predicted,  although  the  magnitudes  are  significantly 
different.  Thus,  the  reacting,  recirculating  flow  has  not  been  satisfactorily  predicted. 
Although  the  chemistry'  model  used  in  the  present  study  is  undeniably  crude,  the 
computational  experience  gained  in  this  study  implies  that  the  principal  effort  in  the 
future  should  be  to  develop  a turbulence  transport  model  that  will  accurately  predict  the 
near  field  of  the  jet  mixing  process  in  this  class  of  recirculating  flows.  The  calculation 
time  for  this  flow  field  was  about  40  min. 

5.0  CONCLUDING  REMARKS 

The  calculation  procedure  for  recirculating  combustor  flows,  based  on  solving  a 
system  of  equations  for  vorticity.  stream  function,  and  species  parameter  by  using  the 
decay-function,  finite-difference  scheme  of  Chien,  has  been  evaluated  by  comparisons 
with  four  sets  of  experimental  data  - turbulent  pipe  flow,  and  three  recirculating  flows. 
The  turbulent  pipe  flow  (Laufer)  and  the  case  of  constant  density'  recirculating  flow 
(Barchilon  and  Curtet)  were  predicted  with  acceptable  accuracy.  The  axial  velocity  field 
of  the  nonreacting,  variable-density  recirculating  flow  (Schulz)  was  acceptably  well 
predicted,  but  the  predicted  species  mass  fraction  field  showed  too  great  of  a rate  of 
diffusion  and  spreading.  In  this  calculation,  the  species  were  conserved  to  within  eight 
percent  of  the  experimental  results  in  the  far  field.  Neither  the  axial  velocity  field  nor 
the  species  field  of  the  reacting,  variable-density  combustor  flow  (Chriss)  were  well 
predicted.  The  theoretical  calculations  showed  a large  velocity  overshoot,  but  it  could  not 
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be  determined  whether  the  viscosity  model  or  the  chemistry  model  caused  this  effect. 
The  far  field  species  concentration  was  75  percent  greater  than  experiment,  indicating 
that,  for  this  case,  the  solution  did  not  conserve  species. 

It  was  shown  that,  for  variable-density  recirculating  flows,  the  vorticity  source  term 
should  be  complete  because  the  terms  that  have  been  considered  negligible  in  other 
studies  have  a significant  effect  on  the  quality  of  the  solution  in  the  near  field.  Also,  it 
was  found  that  because  of  the  retained  source  terms,  the  inlet  boundary  values  of 
velocity,  density,  and  viscosity  must  be  consistent  because  they  affect  the  numerical 
solution  downstream  by  interaction  with  the  source  terms  in  the  vorticity  equations. 

The  primary  research  goal  in  the  future  studies  should  be  to  develop  turbulent 
viscosity  models  that  are  applicable  to  the  strong,  developing  shear  layers  that  exist  in  the 
recirculating  flows  of  the  type  studied  herein.  Until  an  adequate  viscosity  model  is 
obtained,  it  is  probably  not  worthwhile  considering  more  complex  chemistry  models  for 
these  flows. 


REFERENCES 

1.  Schulz,  R.  J.  "An  Investigation  of  Ducted,  Two-Stream,  Variable  Density,  Turbulent 

Jet  Mixing  with  Recirculation."  AEDC-TR-76-152  (ADA034537),  AFOSR-TR- 
76-1087,  January  1977. 

2.  Schulz,  R.  J.  "An  Experimental  and  Theoretical  Investigation  of  Confined, 

Two-Stream  Variable  Density,  Turbulent  Jet  Mixing  with  Recirculation."  Ph.D. 
Dissertation,  The  University  of  Tennessee,  Knoxville,  1976. 

3.  Gosman,  A.  D.,  Pun,  W.  K.,  Ruchal,  A.  K.,  Spalding,  B.  D.,  and  Wolfshtein,  M.  Heat 

and  Mass  Transfer  in  Recirculating  Flows.  Academic  Press,  New  York,  1969. 

4.  Chien,  J.  C.  "Numerical  Analysis  of  Turbulent  Separated  Subsonic  Diffuser  Flows." 

AEDC-TR-76-159  (ADA036005),  February  1977. 

5.  Bird,  R.  B..  Stewart,  W.  E.,  and  Lightfoot,  E.  N.  Transport  Phenomena.  John  Wiley 

and  Sons,  Inc.,  New  York,  1960. 

6.  Schlichting,  H.  Boundary  Layer  Theory'.  Sixth  Edition,  McGraw-Hill  Book  Company, 

Inc.,  New  York,  1966. 

7.  Launder,  B.  E.,  and  Spalding,  D.  B.  Lectures  on  Mathematical  Models  of  Turbulence. 

Academic  Press,  Inc.,  New  York,  1972. 


23 


A EDC-TR-78-29 


8.  Elghobashi,  S.  E.  "Characteristics  of  Gaseous  Turbulent  Diffusion  Flames  in 

Cylindrical  Chambers."  Ph.D.  Dissertation,  Mechanical  Engineering  Department, 
Imperial  College  of  Science  and  Technology.  London,  England,  1974. 

9.  Chien,  J.  C.  "A  General  Finite-Difference  Formulation  with  Application  to 

Navicr-Stokes  Equations."  Journal  of  Computational  Physics,  Vol.  20,  Number 
3,  March  1976,  pp.  268-278. 

10.  Chien,  J.  C.  "A  General  Finite-Difference  Formulation  with  Application  to 

Navicr-Stokes  Equations."  Computers  and  Fluids,  Vol.  5,  1977,  pp.  15-31. 

11.  Laufer.  J.  "The  Structure  of  Turbulence  in  Fully  Developed  Pipe  Flow."  NACA 

Report  1174,  1954. 

12.  Barchilon,  M.  and  Curtet,  R.  "Some  Details  of  the  Structure  of  an  Axisymmetric 

Confined  Jet  with  Backflow."  Transactions  of  the  ASME  Journal  of  Basic 
Engineering,  Vol.  86,  December  1964,  pp.  777-787. 

13.  Becker,  H.  A.,  Mottle,  H.  C..  and  Williams.  G.  C.  "Mixing  and  Flow  in  Ducted 

Turbulent  Jets."  Proceedings  of  the  Ninth  Symposium  (International)  on 
Combustion,  Academic  Press.  Inc.,  Vol.  I,  New  York,  1963,  pp.,  7-20.. 

14.  Chriss,  D.  E.  "An  Experimental  Investigation  of  Ducted,  Reactive,  Turbulent  Jet 

Mixing  with  Recirculation."  AEDC-TR-77-56  (ADA044110),  AFOSR-TR- 
77-0749,  September  1977. 


24 


Radial 

Direction 


to 

U\ 


us  is  the  secondary  jet  velocity  at  nozzle  exit. 

Mean  streamline  pattern  shown  in  upper  half  figure. 
Mean  axial  velocity  profiles  shown  in  lower  half  figure. 

Figure  1.  Turbulent,  ducted,  mixing  system  with  recirculation. 
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Figure  2.  Velocity  profile  across  central  jet  nozzle  used 
to  define  u,  fi,  and  i//  distributions  for  the 
ducted  jet-mixing  calculations. 


Experiment,  Laufer  (Ref.  11) 

o Prandtl  Mixing  Length  Viscosity  Model;  Cp  = 0. 09 


Figure  3.  Comparison  of  experimental  and  theoretical  distributions 
of  axial  velocity  and  turbulent  kinetic  energy  across  a 
turbulent  pipe  flow. 
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Experiment,  Laufer  (Ref.  11) 

o Prandtl  Mixing  Length  Viscosity  Model;  Cp  =0.09 
Law-of-the-Wall  Matching  tor  Boundary  Conditions 


Figure  4.  Comparison  of  experimental,  and  theoretical  distributions 
of  axial  velocity  and  turbulent  kinetic  energy  across  a 
turbulent  pipe  flow  with  law-of-the-wall  matching. 
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Law-of-the-Wall  Matching  for  Boundary  Conditions 


Figure  5.  Comparison  of  experimental  and  theoretical  distributions 
of. axial  velocity  and  turbulent  kinetic  energy  across  a 
turbulent  pipe  flow  with  K-e.  viscosity  model. 
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Experiment,  BarchilonandCurtet  (Ref.  12),  Cy*  0.305 

o Theory 

□ Theory  With  S p (Eq.  (47))  Deleted  from  Eq.  (3) 


Axial  Distance,  Z/Rp 


Figure  6.  Comparison  of  experimental  and  theoretical  distributions 
of  axial  velocity  on  the  centerline  for  constant-density, 
ducted  jet  mixing. 
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Figure  7.  Comparison  of  experimental  and  theoretical  velocity 
distributions  at  three  axial  locations  for  constant- 
density,  ducted  jet  mixing. 
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Axis  Wall 

Experiment,  Schulz  (Ref.  1) 

o □ Theory 

• Theory  with  S0  (Eq.  (47))  Deleted  from  Eq.  (3) 
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Figure  8.  Comparison  of  experimental  and  theoretical  distributions 
of  axial  velocity  and  hydrogen  mass  fraction  for 
variable-density,  ducted  jet  mixing. 


Experiment,  Schulz  (Ref.  1) 

o Theory 
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Figure  9.  Comparison  of  experimental  and  theoretical  velocity 
distributions  at  three  axial  locations  for 
variable-density,  ducted  jet  mixing. 
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Experiment,  Schulz  (Ret.  1) 

o Theory 


Experiment  Theory 
Z/Rq  0.38  0.46 


Experiment 
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Experiment  Theory 
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Figure  10.  Comparison  of  experimental  and  theoretical  distributions 
of  hydrogen  mass  fraction  at  three  axial  locations  for 
variable-density,  ducted  jet  mixing. 
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Experiment,  Chriss  (Ref.  14) 

o a Theory 

(3)  Velocity  Measured  by  Laser  Velocimetry 

(D  Velocity  Calculated  from  Measured  Total  and  Static  Pressures  and  Hydrogen  Atom 
Mass  Fractions  and  Based  on  An  Equilibrium  Chemistry  Hypothesis 


Figure  11.  Comparison  of  experimental  and  theoretical  distributions 
of  axial  velocity  and  hydrogen  mass  fractions  for 
variable-density,  reacting  jet  mixing. 
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Figure  12.  Comparison  of  experimental  and  theoretical  velocity 
distributions  at  three  axial  locations  for 
variable-density,  reacting  jet  mixing. 
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Figure  13.  Comparisons  of  experimental  and  theoretical  distributions 
of  hydrogen  atom  mass  fractions  at  three  axial  locations 
for  variable-density,  reacting  jet  mixing. 
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APPENDIX  A 

COEFFICIENTS  IN  THE  STANDARD  FORM  OF  THE  GOVERNING  EQUATION 
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IV.  Turbulent  Kinetic  Energy 


V.  Dissipation 
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NOMENCLATURE 


sl 0,  h2^ 
b2^ 

C 

**** 

c 


Cp 

CT 

d0,  D0 

G0 

Href 

I 

J 

K 

e 


ni|. , m<p,  mp,  mi 
P 


r 

r 


Rd 

Rg 

Ri 

T 


Coefficients  in  the  standard  form  of  the  governing  equations,  Eq.  (38) 
Coefficients  in  the  standard  form  of  the  governing  equations,  Eq.  (38) 
Species  parameter,  either  mass  fraction  of  air,  or  given  by  Eq.  (22) 
Species  parameter,  Eq.  (23) 

Specific  heat  at  constant  pressure 

Craya-Curtet  Number,  Ref.  1 

Source  term  for  the  dependent  variable 

Decay  function  for  the  dependent  variable  0,  Eq.  40 

Reference  enthalpy 

Mesh  index  for  axial  position 

Mesh  index  for  radial  position 

Turbulent  kinetic  energy  per  unit  mass 

Prandtl  Mixing  Length,  Eq.  (11) 

Mass  fractions  of  fuel,  oxidizer,  product  and  inert  species,  respectively 
Static  pressure 
Radial  coordinate 

Nondimensional  radial  coordinate,  (Rp  - r)/Rp , zero  at  the  wall 

Duct  radius 

Universal  gas  constant 

Jet  radius 

Static  temperature 
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U 

U 

Uj 

Um 

US 

Uw 

v* 

v 

W 

z 

6* 

e 

f 

T? 

M , 
Mt 

V 

P 

* 

<t> 

Tw 


Axial  velocity 

Nondimensional  axial  velocity,  u = u/um 
Jet  velocity 

Maximum  axial  velocity  at  a given  axial  position 

Secondary  stream  velocity 

Wall  velocity  at  7 = 0.05 

Friction  velocity,  V*  = ^rw/p 

Radial  velocity 

Molecular  weight,  Eqs,  (18)  or  (26) 

Axial  coordinate 

Boundary  layer  displacement  thickness 

Dissipation  rate  of  turbulent  kinetic  energy 

Stretched  axial  coordinate 

Stretched  radial  coordinate 

Dynamic  viscosity,  ti  = pt  + juL 

Turbulent  viscosity 

Laminar  viscosity 

Kinematic  viscosity,  v = Ip 

Density 

Stream  function,  Eq.  (5) 

Dependent  variable,  either  \jj.  C,  K,  or  e 
Wall  shear  stress 
Vorticity,  Eq.  (3) 
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SUBSCRIPT 

FS 

JN 

RS 

* 

V 

SUPERSCRIPT 


Forward  stagnation  point 

Radial  matching  point  outside  the  laminar  sublayer  in  pipe  flow 

Rear  stagnation  point 

With  respect  to  stretched  axial  coordinate 

With  respect  to  stretched  radial  coordinate 


An  over  bar  represents  a time-averaged  quantity 
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